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1. Introduction 

oo 

' The Green's function method is a very convenient formalism in condensed mat- 

ter physics, and many progresses have been achieved in the last fifty years. When 
applied to interacting systems, such an approach is usually based on the hypothesis 
that the interaction among the particles is weak and can be treated in the frame- 
work of some perturbation schemes. In this line of thinking a consolidated scheme 
has been constructed, mostly based on diagrammatic expansions, Wick's theorem, 
Dyson equation, and so on. However, in the last few decades new materials with 
unconventional properties have been discovered. It is believed that the origin of 
such anomalous behaviors is generally due to strong electronic correlations in nar- 
row conduction bands [1]. In this line of thinking many analytical methods have 
been developed for the study of strongly correlated electron systems [2]. The main 
difficulties are connected with the absence of any obvious small parameter in the 
strong coupling regime and with the simultaneous presence of itinerant and atomic 
aspects. The concept that breaks down is the existence of the electrons as particles 
with some well-defined and intrinsic properties. The presence of interaction modifies 
the properties of the particles: what are observed are new particles with new pecu- 
liar properties entirely determined by the dynamics and by the boundary conditions. 
These new objects appear as the final result of the modifications imposed by the 
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interactions on the original particles and contain, by the very beginning, the effects 
of correlations. The choice of new fundamental particles, whose properties have to 
be self-consistently determined by dynamics, symmetries and boundary conditions, 
becomes relevant. 

As a simple example, let us consider an atomic system described by the Hamil- 
tonian 

h = -h ^2 riv* + V( fWWm ( L1 ) 

cr 

(Pa- denotes an Heisenberg electronic field with spin a =f, |, satisfying canonical 
anticommutation relations; \i is the chemical potential and V is the strength of the 
interaction. This model is exactly solvable in terms of the operators 

£ CT = ip a ip_ a ipl a T] a = ip a ipl a ip_ a (1.2) 

which are eigenoperators of the Hamiltonian 

i-C = K, H] = -rt \-n = [77, H] = -Qi - V)r, (1.3) 

Due to the presence of the interaction, the original electrons ip a are no more 
observables and new stable elementary excitations, described by the field operators 
£ and rj, appear. Due to the ^-interaction, two sharp features develop in the band 
structure: the energy level E = —fi of the bare electron splits in the two levels 
Ei = —pi and E 2 = V — \x. The bare electron reveals itself to be precisely the wrong 
place to start. A perturbative solution will never give the band splitting. 

On the basis of this evidence one can be induced to move the attention from 
the original fields to the new fields generated by the interaction. The operators 
describing these excitations, once they have been found, can be written in terms of 
the original ones and are known as composite operators. 

The convenience of developing a formulation to treat composite excitations as 
fundamental objects has been noticed for the many-body problem of condensed 
matter physics since long time. Recent years have seen remarkable developments in 
many-body theory in the form of an assortment of techniques that may be termed 
composite particle methods. The beginnings of these types of techniques may be 
traced back to the work of Bogolubov [3] and later to that of Dancoff [4]. The work 
of Zwanzig [5], Mori [6] and Umezawa [7] has to be mentioned. Closely related to 
this work is that of Hubbard [8-10], Rowe [11], Roth [12] and Tserkovnikov [13, 14]. 
The slave boson method [15-17], the spectral density approach [18,19] and the 
composite operator method (COM) [20-27] are also along similar lines. This large 
class of theories is founded on the conviction that an analysis in terms of elementary 
fields might be inadequate for a system dominated by strong interactions. 

All these approaches are very promising because all the different approximation 
schemes are constructed on the basis of interacting particles: some amount of the in- 
teraction is already present in the chosen basis and permits to overcome the problem 
of finding an appropriate expansion parameter. However, one price must be paid. 
In general, the composite fields are neither Fermi nor Bose operators, since they 
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do not satisfy canonical (anti) commutation relations, and their properties, because 
of the inherent definition, must be self-consistently determined. They can only be 
recognized as fermionic or bosonic operators according to the number, odd or even, 
of the constituting original electronic fields. New techniques of calculus have to be 
developed in order to treat with composite fields. In developing perturbation cal- 
culations, where the building blocks are now the propagators of composite fields, 
the consolidated scheme (diagrammatic expansions, Wick's theorem, . . . ) gives rise 
to very complicated approaches [28] whose application is far from being easy. The 
formulation of the Green's function method must be revisited and new frameworks 
of calculations have to be formulated. 



2. Green's Function and Equation of Motion Formalism 

Let us consider a system of N e interacting Wannier-electrons residing on a Bra- 
vais lattice of iV sites, spanned by the vectors Ri = i. We ignore the presence of 
magnetic impurities and restrict the analysis to single-band electron models. The 
generalization of the formalism to more complex systems is straightforward [see for 
example [29-33]]. The system is enclosed in a finite but macroscopically large volume 
V and is supposed to be in a thermodynamical equilibrium state at a finite tem- 
perature T. In a second quantization scheme this system is described by a certain 
Hamiltonian 

H = H[tp(i)] (2.1) 

describing, in complete generality, the free propagation of the electrons and all the 
interactions among them and with external fields. ip(i) denotes an Heisenberg elec- 
tronic field [ i — (i, t)] satisfying canonical anticommutation relations derived from 
the Pauli principle 

{(Pa{ht), ( pl(},t)} = S iJ S <Ta K(M),</Uj,0} = (2.2) 

Any physical property of the system can be connected to the expectation value 
of a specific operator. The expectation value of an arbitrary operator A = A[ip(i)] 
can be computed, for the grand canonical ensemble, by means of 

(A) = ^[e-**-^] 3 

V 1 y r [ e -/3(ff- M JV)] V • ; 

where the trace implies a sum over a complete set of states in the Hilbert space. 
N = J2i a V^J-WwW is the total number operator, (3 = (fc^T) -1 is the inverse 
temperature, \i is the chemical potential which is fixed in order to give the chosen 
average number of particles N e = (N). 

To evaluate the expectation value (A), it is possible to use the equation of motion 

i^) = W),H] (2.4) 
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in order to derive one or more equations for this quantity or, better, for the corre- 
sponding Green's functions, as explained below. However, the equation of motion 
(2.4) generates higher-order operators and more and more complex equations are 
needed. The traditional approximation schemes, often based on perturbative cal- 
culations, use as building blocks the noninteracting propagators. The mean-field 
formulation, which corresponds to a linearization of the equation of motion (2.4), 
also belongs to this category. 

On the hypothesis that the original fields are not a good basis, we choose a set of 
composite fields {ip{i)} in terms of which a perturbation scheme will be constructed. 
Firstly, we choose the set according to the physical properties we want to 
study. Roughly, the properties of electronic systems can be classified in two large 
classes: single particle properties, described in terms of fermionic propagators, and 
response functions, described in terms of bosonic propagators. These two sectors, 
fermionic and bosonic, are not independent but interplay each other, and a fully 
self-consistent solution usually requires that both sectors are simultaneously solved. 
Once the sector, fermionic or bosonic, has been fixed, we have several criteria for 
the choice of the new basis. In constructing the composite fields no recipe can be 
given without thinking to its drawbacks, but many recipes can assure a correct 
and controlled description of relevant aspects of the dynamics. One can choose: the 
higher order fields emerging from the equations of motion (i.e., the conservation of 
some spectral moments is assured), the eigenoperators of some relevant interacting 
terms (i.e., the relevant interactions are treated exactly), the eigenoperators of the 
problem reduced to a small cluster, . . . 

Let be a n-component field 



(2.5) 



We do not specify the nature, fermionic or bosonic, of the set In the case of 

fermionic operators it is intended that we use the spinorial representation 

^ = ( fc« ) = ^ vrM ^ m(<) ) (2 - 6) 

The dynamics of these operators is governed by the given Hamiltonian H = 
H[ip(i)] and can be written as 

i^(<) = [m,H] = J{i) (2.7) 

It is always possible to decompose the source J(i) under the form 

J(i) = e(-iV)ip(i) + 8J(i) (2.8) 

where the linear term represents the projection of the source on the basis {ip} and 
is calculated by means of the equation 

([5j(M),^aow = o (2.9) 
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Here r] = ±1; usually, it is convenient to take rj = 1 (rj = — 1 ) for a fermionic 
(bosonic) set (i.e., for a composite field constituted of an odd (even) number of 
original fields) in order to exploit the canonical anticommutation relations of 
but, in principle, both choices are possible. Accordingly, we define 

_( {A,B} = AB + BA for V = 1 
LA ^ ~ \ [A, B] = AB — BA for rj = -1 1 UJ 

(■■•) denotes the quantum statistical average over the grand canonical ensemble, 
according to Eq. (2.3). Hereafter, unless otherwise specified, time and space trans- 
lational invariance will be considered. The action of the derivative operator e(— iV) 
on is defined in momentum space 

e(-iV)iKi) =e(-iV)-^^e lk ' R ^(k,t) = -^^e lk ^(k)^(k,t) (2.11) 

where k runs over the first Brillouin zone. The constraint (2.9) gives 

m(k) = e(k)/(k) (2.12) 
after defining the normalization matrix 

J(i,j) = ([^(i,t),^(j,t)U = ^e ik ( R '-^7(k) (2.13) 

k 

and the m-matrix 

m(i, j) = ([J(i, t), ^(j, 01,> = ^ E e ik ' (Ri ~ Rj) m(k) (2.14) 

k 

Since the components of ip{i) contain composite operators, the normalization 
matrix J(k) is not the identity matrix and defines the spectral content of the ex- 
citations In fact, the use of composite operators has the advantage of describing 
crossover phenomena as the phenomena in which the weight of some operator is 
shifted to another one. 

It is worth noting that the normalization matrix 7(k) and the m(k) matrix are 
the lowest order generalized spectral moments M^(k) which are defined as 

M«(k) = F.r.([(ia/9t)^(i J o,^ t a*)] fJ ) (2-15) 

where FT. stands for the Fourier transform. The generalized spectral moments 
M^ p \k) have the Hermiticity property 

M«(k) = (M£>(k))* (2.16) 

since at the equilibrium M^ik) is time-independent. 
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Coming back to our original problem, the evaluation of the expectation value (A), 
it is possible to use the equation of motion ij^(i) = [ijj(i),H] in order to derive 
equations for (A). However, the correlation functions satisfy homogeneous equations. 
A convenient generalization of the concept of correlation functions is furnished by 
the Green's functions (GF) which have some advantages in the construction and 
solution of the equations determining them. In particular, the two-time Green's 
functions contain most of the relevant information on the properties of the system: 
expectation values of the observables, excitation spectrum, response to external 
perturbation, and so on. Different types of GF can be defined; for statistical systems 
it is better to consider the real time thermodynamic Green's functions where the 
averaging process of the Heisenberg operators is performed over the grand canonical 
ensemble. 

By considering the two-time thermodynamic Green's functions [34-36], let us 
define the causal function 



G c (i, 3 ) = (CWWti)]) =9(t l -t 3 ){^i)^U)) -rfl(tj -tiX^OX*)) (2-17) 



G R ' A (i,j) = {R,A[1>(i)tfti)]) = ±e[±{t l -tMm)^m,) (2.18) 



By means of the Heisenberg equation (2.7) and using the decomposition (2.8), 
the Green's function G Q (i,j) = (Q[ip(i)^(J)]), where Q = C, R, A, satisfies the 



A($)G*(i, j^idj) = A(9 i )G?(i,i)At(9 j ) + {Q[5J{i)5J\j)\) (2.19) 



the retarded and advanced functions 



equation 



where the derivative operator A(<9j) is defined as 



m) 



.d_ 



e(-iV0 



(2.20) 



and the propagator Gq (i, j) is defined by the equation 

HdJG$(i,j) = i8(t i -t j )I{i,i) 



(2.21) 



By introducing the Fourier transform 




(2.22) 



equation (2.19) in momentum space is written as 




(2.23) 



where the self-energy S^*(k, uo) has the expression 

£ Q *(k,cu) = I- 1 (k)B Q (k,to)r 1 (k) 



(2.24) 
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with 

B Q {k,co) = F.T.{Q[5J{i)5J\j)\) (2.25) 
Next, we introduce the irreducible self-energy E^(k, a;) by means of the definition 

£ Q (k, co)G Q ik, co) = J(k)£ Q *(k, w)G#(k, co) (2.26) 

Equation (2.23) takes the form 

G Q {k, co) = G®(k, co) + C7^(k, cj)/- 1 (k)S Q (k, co)G Q (k, co) (2.27) 

and can be formally solved as 

G Q (k, co) = — 1 Z(k) (2.28) 

co — e(k) — L^(k, co) 

The formal definition (2.24) of self-energy must be manipulated to avoid any 
flowing on tautology. By noting that C7^ _1 (k, co) = / _1 (k)[u; — e(k) — X^(k, co)] we 
can express £ Q (k, co) as 

YP{k,u) = B% r {k,u)r\k) (2.29) 
where B2 r (k, co) indicates the irreducible part of the propagator B®(k, co), given by 

B2 r (k,co) = (2.30) 

irrV ' fiQ- 1 (k,cu) + /- 1 (k)C7j(k,cu)/- 1 (k) 

We have constructed a generalized perturbative approach designed for formu- 
lations using composite fields. Equation (2.27) is a Dyson-like equation and may 
represents the starting point for a perturbative calculation in terms of the propaga- 
tor C7^(k, co). Contrarily to the usual perturbation schemes, the calculation of the 
"free propagator" C7^(k, co) is not an easy task and large part of this article will 
be dedicated to this problem. Then, the attention will be given to the calculation 
of the self-energy S^(k, co), and some approximate methods will be presented. It 
should be noted that the computation of the two quantities G® (k, co) and S Q (k, co) 
are intimately related. The total weight of the self-energy corrections is bounded by 
the weight of the residual source operator 5J(i). According to this, it can be made 
smaller and smaller by increasing the components of the basis ip(i) [e.g. by including 
higher-order composite operators appearing in 8J(i)\. The result of such a procedure 
will be the inclusion in the energy matrix of part of the self-energy as an expansion 
in terms of coupling constants multiplied by the weights of the newly includes basis 
operators. In general, the enlargement of the basis leads to a new self-energy with a 
smaller total weight. However, it is necessary pointing out that this process can be 
quite cumbersome and the inclusion of fully momentum and frequency dependent 
self-energy corrections can be necessary to effectively take into account low-energy 
and virtual processes. According to this, one can chose a reasonable number of com- 
ponents for the basic set and then use another approximation method to evaluate 
the residual dynamical corrections. 



7 



F. Mancini and A. Avella 



3. GF properties, spectral representation, zero-frequency func- 
tions 

3.1. Equations of motion 

In the previous Section we have constructed a generalized perturbative approach 
based on a Dyson equation designed for formulations using composite fields. Two 
quantities appear in the Dyson equation (2.27): the "free propagator" G^(k, u) and 
the self-energy S^(k, uS). By postponing to next Sections the problem of computing 
the self-energy, in this Section we concentrate on the calculation of the Green's 
functions Go (k, cu) which constitute the building blocks of the perturbation scheme 
we are trying to formulate. For the sake of simplicity we will drop the sub index 
in the definition of G®(k, u). 

One fundamental aspect in a Green's function formulation is the choice of the 
representation. The knowledge of the Hamiltonian and of the operatorial algebra is 
not sufficient to completely specify the GF. The GF refer to a specific representation 
(i.e., to a specific choice of the Hilbert space) and this information must be supplied 
as a boundary condition to the equations of motion that alone are not sufficient to 
completely determine the GF. As well known, the same system can exist in different 
phases according to the external conditions; the existence of infinite inequivalent 
representations [37-39] where the equations of motions can be realized, allows us 
to pick up, among the many possible choices, the right Hilbert space appropriate 
to the physical situation under study. The use of composite operators leads to an 
enlargement of the Hilbert space by the inclusion of some unphysical states. As a 
consequence of this, it is difficult to satisfy a priori all the sum rules and, in general, 
the symmetry properties enjoined by the system under study. In addition, since the 
representation where the operators are realized has to be dynamically determined, 
the method clearly requires a process of self-consistency. 

From this discussion it is clear that fixing the representation is not an easy 
task and requires special attention. In the literature the properties of the GF are 
usually determined by starting from the knowledge of the representation. Owing to 
the difficulties above discussed we cannot proceed in this way. Therefore, we will 
derive the general properties of the GF on the basis of the two elements we have: 
the dynamics, fixed by the choice of the Hamiltonian (2.1), and the algebra, fixed 
by the choice of the basic set (2.5). The problem of fixing the representation will be 
considered in the next Sections. 

Let ip(i) be a n-component field satisfying linear equations of motion 



with the energy matrix j) defined by (2.12)-(2.14). If the fields ifj(i) are eigenop- 
erators of the total Hamiltonian, the equations of motion (3.1) are exact. There are 
many non-trivial realistic systems for which it is possible to obtain a complete set of 
eigenoperators of the Hamiltonian (for instance see [40-43]). If the fields ip(i) are not 




(3.1) 
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eigenoperators of H, the equations are approximated; they correspond to neglecting 
the residual source operator 6J(i) in the full equation of motion [see (2.7) and (2.8)] 
and all the formalism is developed with the intention of using the propagators of 
these fields as a basis to set up a perturbative scheme of calculations on the ground 
of the Dyson equation (2.27) derived in the previous Section. 

By means of the field equation (3.1) the Fourier transforms of the various Green's 
functions defined by (2.16), (2.17) and (2.18) satisfy the following equation 

[uj-e(k)]G QM (k,uj) =/(">(k) (3.2) 

where the dependence on the parameter rj has been explicitly introduced. As men- 
tioned in Section 2, the set can be fermionic or bosonic and the parameter r] 
generally takes the value 77 = 1 (77 = —1) for a fermionic (bosonic) set The 
three Green's functions G c , G R and G A satisfy the same equation of motion which 
alone is not sufficient and must be supplemented by other equations. Indeed, the 
GF are determined by solving a first order differential equation of motion, thereby 
the GF are given only within an arbitrary constant of integration. The retarded and 
advanced GF can be completely determined because the factor 9[±(ti — tj)] provides 
the boundary condition: G R,A (i,j) = for t; L = tj =p 5. The determination of the 
causal GF is not so immediate. The most general solution of equation (3.2) is 

C7^)(k,c) = E{p(^^) -i^[o;-^(k)]^(k)} (3.3) 
where c^(k) are the eigenvalues of the e(k), cr^(k) are defined by 

n 

a«(k) = V al (k)n; c \k)I cb (k) a,b=l,...,n (3.4) 

c=l 

f2(k) is the nxn matrix, whose columns are the eigenvectors of e(k). We note that 
while the spectral density matrix erW(k) is completely determined by the energy 
e(k) and normalization 7^(k) matrices, the matrix g^'^ik) is not fixed by the 
equations of motion and must be determined by means of the boundary conditions. 
P represents the principal value. 

By recalling the retarded and advanced nature of G R ' A ^\i,j) it is immediate to 
see that 

g ^>) R (k) = -g^ A (k) = <7<''">(k) (3.5) 

Then, the retarded and advanced functions are completely determined in terms of 
the matrices e(k) and J^(k). As well known, as functions of cu the G R,A ^ (k, uj) are 
analytic in the upper and lower half-planes, respectively. 

The determination of g^ l,r ^ c {k) requires more work. From the definitions (2.16), 
(2.17) and (2.18) we can derive the following exact relations 

G R(n) {iJ) + G Mn) {iJ) = 2G c ^(i,j) - <[^(i), (3.6) 
G R W(iJ)-G A W(iJ) = (3-7) 
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where there appear the two correlation functions 

c^(ij) = (m^u)) c^(ij) = (^ma)) (3.8) 

The definition of C^^(i,j) must be interpreted at level of matrix elements: 
C^^. ha {i,j) = (ipl(j)ijj a (i)}; no scalar product is intended, neither in the spin space. 
These functions are defined for all real times and the equations of motion they obey 
contain no inhomogeneous terms involving a delta function. Indeed, by means of the 
equations of motion (3.1) the Fourier transform of these correlation functions satisfy 
the homogeneous equations 

[u - £(k)]C^ t (k, u) = [to- s(k)]C^(k, u) = (3.9) 

These equations tell us that the Fourier transforms of the correlation functions are 
zero unless the frequency u> is equal to one of the energy levels ui(k) of the system. 
The solutions of (3.9) have the general form 

n 

C^(k,u) = £*[w-w,(k)] C » t (k) (3.10) 
i=i 

n 

C^(k,u) = £*[ W -w,(k)] C 2 (k) (3.11) 



i=i 



with the momentum-dependent Fourier components (k) and S- 1 ] (k) to be de- 
termined. 

We now recall the Kubo-Martin-Schwinger (KMS) relation 

(A(t)B(t)) = {B(t)A(t + i(3)) (3.12) 

where A(t) and B(t) are Heisenberg operators at time t. This relation implies that 
the Fourier transforms C^t(k, uj) and C^ty,(k, u>) are related 

C^(k,u) = e~ puJ C^{k,uj) (3.13) 

and the 77— commutator ([tp(i),^(j)]n) can be expressed in terms of the correlation 
function as 

mUHm = / ^e ik ^)- i ^-^[l+7 7 e-^]C^t(k,a;) (3.14) 

By putting (3.14) into (3.6) and (3.7) and by taking into account Eqs. (3.3), 
(3.5) and (3.10) we obtain 

- W| (k)] L^ c (k) - Ul - rje-^cfjk)} = (3.15) 
1=1 ^ > 

jS[o; - u;,(k)] |^' T/) (k) - hi + ?7 e-^] C « t (k)} = (3.16) 
1=1 ^ J 

The solution of Eqs. (3.15) and (3.16) is remarkably different according to the 
value of the parameter 77 and we shall treat separately the two cases. 
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3.2. Fermionic fields 

For the case of fermionic fields it is convenient to choose 77 = 1. Then, the solution 
of Eqs. (3.15) and (3.16) is 



C (0 (k) =7T 



1 + tanh 



/Mk) 



a^ +1 )(k) 



g^ c \k) = tanh (^^j * {l > +1 \k) 



(3.17) 
(3.18) 



By putting (3.5) and (3.17) into (3.3), (3.10) and (3.11) we get the following 
general expressions for the Green's functions and correlation functions 



C7^ +1 )(k,c) = ^-f 

n 



(k) 



1=1 

n 



UJ - UJlik) ± l8 

l-f F [u t (k)] / F k(k)] 



n 

C^(k, w) = 7T 5 [a; - wj(k)] 



a; — a>j(k) + i5 — a>j(k) — 18 
/3wj(k 



1 + tanh 



1 — tanh 



^(k) ^j 



(7^ +1 )(k) 
(i,+l)( k ) 



(3.19) 

(3.20) 
(3.21) 
(3.22) 



where /f(<^) is the Fermi distribution function: /f(<^) = e ^ +1 ■ 

By recalling that all the fermionic energies are shifted by the chemical potential, 
the locus in the /c-space, defined by uJi(k) = 0, will define the Fermi surface. By 
looking at (3.15b) we can see that the imaginary part of the causal GF vanishes on 
the Fermi surface. In the fermionic case the right procedure of calculation is to start 
from the retarded (advanced) GF, and then to compute the other Green's functions 
and correlation functions by means of the relations 



&[G c( - +1 \k,u)] = $t[G R ' A ( +1 \k,u;)} 



3[C7 c ( +1 )(k,cj)] = ±tanh 



C^t (k, u) = =F 



1 + tanh 



(3uj 

IT 



%[G R > A ( +1 \k,u)} 
Q[G R ^ +1 \k,uj)] 



We note the dispersion relations 

1 p+00 I 
?H[G R ' A{+1 \k, u)] = T-P / duo — —^s[G R ^ +1 \k, uj)} 

7T UJ - UJ 



^[G c{+1 \k,iu)] = --P 

71 



-00 
-00 



UJ — UJ 



C oth(^-)%G ci * +l \k,uj)} 



(3.23) 
(3.24) 

(3.25) 

(3.26) 
(3.27) 
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By introducing the spectral function 

n 

p^(k,co) = ^5[c-c,(k)]^ +1 )(k) = T -^[G R ^ +1 \k,u)\ 

* — ' 7T 

1=1 

we can establish the spectral representation 



(3.28) 



„ to — uj' ± id 



/+oo 
dJ p { + x \k,J 
-oo 



1-/p(o/) + / P (u/) 



a; — a;' + iS uj — uj' — \5 



(3.29) 
(3.30) 



3.3. Bosonic fields 



For the case of bosonic fields it is convenient to choose 77 = — 1. For any given 
momentum k we can always write 



wj(k) 



= for/G A(k) CJV = {l,...,n} 
^ for/G £(k) — N — A(k) 



(3.31) 



Obviously, A(k) can also be the empty set (i.e., A(k) = and B(k) = N). For 
/ G B(k) the solution of (3.15) and (3.16) is 



^t( k ) = 7r 



1 + coth 



g V,-i)c {k ) = coth a^ik) V/ G S(k) 



(3.32) 
(3.33) 



For / G -<4(k) the Fourier coefficients (k) and g^ 1, ^ c (k) cannot be determined 
from Eqs. (3.16). It is convenient to introduce the function T(k) 



2,71 

«eA(k) 



(3.34) 



iGA(k) 



By considering that from (3.16) and (3.10) 



lim[l-e-^C^ t (k,^) = 27r%) V ^'^(k) 

ieA(k) 



(3.35) 



we must distinguish two situations: 

EA(k) 1 



1- IfE^( k) ^'" 1) (k)=0, then 



J](J[w-w,(k)]cW t (k) = 27r(J(w)r(k)+27r ^ (J[w-w,(k)] 



iefi(k) 



e /M(k) _ 1 



^(k) 
(3.36) 
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^ < 5[o;-a;,(k)] C ^(k) = 27r ( 5(a;)r(k)+27r E ^-^(k)] ^ _ ^ ^(k) 

(3.37) 

and 

C7^ t (k,0) = C^(k,0) (3.38) 

2. If Xl«GA(k) <T ^'~ 1 ' ) (k) 7^ but finite, then in order to satisfy (3.35) C^(k,uj) 
must have a singularity of the type ^ in the limit u — > 0. In fact 

(1 - = - \fPJ + -)^ = 1-|*- + - P-39) 

Then 

E <^ - <"i( k )]<£ t ^ = + 2tt £ *[u; - u;, Wl ^wTi ^ 11 ( k ) 

(3.40) 

n n 1 

E ^ - ^( k )]<t (k) = 2 ^H r ( k ) + 27r E ^ - ^( k )]^poii ff( ' rl) ( k ) 

(3.41) 



i=i i=i 
and 



C^ t (k,0)-C^(k,0)= E ^ ,_1) (k) (3.42) 

iGA(k) 

It is clear from (3.40) and (3.41) that the situation where J2ieA(k) c^' -1 ^(k) 7^ 
leads to a situation in which for / e -A(k) the Fourier coefficients (k) and (k) 

diverge as [/3u;;(k)] _1 . Since the correlation function in direct space must be finite, 
at finite temperature this is admissible only in the thermodynamic limit and if the 
dispersion relation u>;(k) is such that the divergence in momentum space is inte- 
grable and the corresponding correlation function in real space remains finite. For 
finite systems and for infinite systems where the divergence is not integrable we 
must have Y2i e ACk) a ^'~^ s ) = 0- The calculation of the spectral density matrices 
(j('' -1 )(k) it not a simple dynamical problem, but requires the self-consistent calcu- 
lation of some expectation values, where the boundary condition and the choice of 
the representation play a crucial role. A finite value of X^eA(k) c^' -1 ^) is generally 
related to the presence of long-range order and the previous statement is nothing 
but the Mermin- Wagner theorem [44]. 

Summarizing, by using (3.5) and by putting (3.36) and (3.37) into (3.15), for 
finite systems and T ^ 0, we get the following general expressions for the GF and 
correlations function 

g r ^ (k, c) = E - ' n I - a ( 3 - 43 ) 

j-^ uj - uji(k) ±i5 
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G^-^k,^) = r(k) 

l€B(k) 



uj + \5 uj — iS 
1 + /bH /b(w) 



C^t(k,u;) = 27rr(k)5(a;)+7r ^ 5[w - w,(k)] 

ZeB(k) 

C^(k,w) = 27rr(k)(JH-7r ]T <S[w-a;,(k)] 

iSB(k) 



— o>j(k) + i<5 — o>j(k) — i5 



(3.44) 



1 + coth 



r('-l) 



(k) 



i-cot h r^ 



w 



(3.45) 
(i,-i)( k) 

(3.46) 

ith the condition that XweA(k) ^'"^(k) = 0. /b(w) is the Bose distribution func- 
tion: /b(^) = e 0^,_ 1 - It is possible to have <7^' _1 )(k) ^ /or I E A(k) only for 
infinite systems and if the divergence is integrable. 

The previous formulas show that when zero-energy modes are present a zero- 
frequency singularity appears in the correlation function C^t (k, a;) and in the imag- 
inary part of the causal function G c '^ 1 - ) (k, m). Such singularity does not contribute 
to the retarded and advanced GF. Then, in the bosonic case the right procedure of 
calculation is to start from the causal GF, and compute the other GF by means of 
the relations 

^[G R ^- l \k,u)] = ^[G c ^(k,uj)} 



9f[G^- 1 )(k,a;)] = ±tanh (^j ^^(k, 



1 + tanh 



f3u> 

IT 



3[G c( - 1} (k, 



uj 



UJ 



We note the dispersion relations 

7T > /_ 00 UJ-UJ 

U[G c ^\k,uj)} = --P [ + °°duj -J— tanh ( ^) Q[G c ^\k,u)] 

7T UJ-UJ \ 2 J 

Also in the bosonic case we can introduce a spectral function 



(3.47) 
(3.48) 

(3.49) 

(3.50) 
(3.51) 



p^iKcu) = , y6[u-u l (k)]a<!- r >(k) = T-^[G RA{ - 1] (k, 

L — / n 



UJ 



(3.52) 



i=i 



However, the zero- frequency function T(k) does not contribute to p' 1 - ) (k, a;) and 
a spectral representation can be established only for the retarded (advanced) GF 



G R ' A <-'\k,uj) = J 



oo 



uj — uj 1 ± i8 



(3.53) 
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For the bosonic causal GF a spectral representation exists only when r(k) = 0: 

1 + /b("') /b(w') 



/+oo 
-oo 



uj — uo' + iS uj — uj' — iS 



(3.54) 



3.4. Sum rules and some useful relations 



Coming back to a generic value of T(k), we note that from the definition (3.4) 
the following sum rule can be derived 

/+oo n 
dup<-"\k,u) = ^V^(k) = I M (k) (3.55) 
-°° i=i 

This is a particular case of a general sum rule. From the results (3.21), (3.22), 
(3.45) and (3.46) we obtain 

am^im = E^E ^-^-^m-^m^ (3 . 5 6) 

1=1 k 

Recalling the expression (2.15) of the generalized spectral moment M^'^(k) we 
immediately have 

/+oo n 
dujuj p p^\k,uj) = ^(k)V ( ' fl) (k) = M (M) (k) (3.57) 
■°° i=i 

Some interesting results can be obtained by noting that the correlation function 
C^ii^j) = (0WV^C?)) an d the energy matrix e(i,j) do not depend on 77. As men- 
tioned above, once we have chosen a basic set fermionic or bosonic, it is a 
only a question of convenience to choose?] = 1 or 77 = — 1. Let us consider the case 
of a bosonic set and let us suppose to perform two series of calculations: one 
with rj = — 1 and another with 77 = 1. Then, it is immediate to obtain the following 
relations: 

1. by equating (3.21) and (3.45) we obtain 

r(k) = \ ^' +1) (k) (3.58) 
ieA(k) 

a^\k) = tanh (^^j a^ +1 \k) V / G B(k) (3.59) 

2. by equating (3.25) and (3.49) we obtain 

SfG^-^k, uj)] = ±%[G R ' A{+1 \k, uj)} (3.60) 
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3. from the sum rule (3.55) and by means of (3.59) and (3.4) we obtain 



^(k) = f>nh f^f±) ^(k^k^'Ck) (3.61) 

^ ' c=l 



1=1 



We see that the general structure of the GF is remarkably different according to 
the statistics. For fermionic composite fields (i.e., when it is natural to choose r] — 1) 
all the Green's functions and correlation functions are completely determined. The 
zero-frequency function T(k), defined on the Fermi surface o>j(k) = /i, contributes 
to the spectral function p( +1 )(k, uS) (see 3.28), it is directly related to the spectral 
density functions <7^' +1 )(k) by means of equation (3.58), and its calculation does 
not require more information. Also, it does not contribute to the imaginary part 
of the causal GF. For bosonic composite fields (i.e., when it is natural to choose 
77 = — 1) the retarded and advanced GF are completely determined, but the causal 
GF and the correlation function depend on the zero-frequency function T(k), defined 
on the surface c^(k) = 0. It is now clear that the causal and retarded (advanced) 
GF contain different information and that the right procedure of calculation is con- 
trolled by the statistics. In particular, in the case of bosonic fields (i.e., for 77 = — 1) 
one must start from the causal function and then use (3.34) to compute the other 
GF. On the contrary, for fermionic fields (i.e., for rj = 1) the right procedure for 
computing the correlation function requires first the calculation of the retarded (ad- 
vanced) function and then the use (3.23), (3.24) and (3.25) to compute the other GF. 
Moreover, it is worth noting that T(k) is undetermined within the bosonic sector 
(i.e., i] = —1). It is true that T(k) could be computed by considering an anticom- 
mutating algebra: remaining in the bosonic sector we make the choice rj — 1 and 
T(k) can be calculated by (3.58) or equivalently by means of the following relation 
T(k) = I lim^^o ujG c( - +1 \\i, uj) which can be easily obtained from (3.44). However, 
the calculation of the cr <7,+1 - ) (k) requires the calculation of the normalization matrix 
/( +1 )(k) that, for bosonic fields, generates unknown momentum dependent correla- 
tion functions whose determination can be very cumbersome as requires, at least in 
principle, the self-consistent solution of the integral equations connecting them to 
the corresponding Green's functions. In practice, also for simple, but anyway com- 
posite, bosonic fields the T(k) remains undetermined and other methods should be 
used. Similar methods, like the use of the relaxation function [45], would lead to the 
same problem. 

Actually, all issues related to T(k) have a natural playground in dealing with 
the ergodicity of the dynamics under investigation. More detail on this topic can be 
found in a manuscript in this same issue [46]. 

The formulation given in this Section needs some modifications in the case of zero 
temperature. In particular, Eqs. (3.15) and (3.16) are not applicable and we must 
proceed in a different way. After a straightforward derivation [26], it is immediate 
to see that the limit T — > of the expressions (3.21), (3.22), (3.45) and (3.46) gives 
the right result. 
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4. A self-consistent scheme 

As stressed in Section 1, in the study of highly interacting systems, where tra- 
ditional perturbative calculations in terms of the noninteracting fields fail, a way 
to reconcile the powerful perturbation theory with the presence of complex and/or 
strong interactions is to describe the system in terms of a new set of fields, com- 
posite operators, generated by the interactions themselves. These field operators 
undoubtedly constitute a better starting point: they appear as the final result of the 
modifications imposed by the interactions on the original particles and contain, by 
the very beginning, the effects of the correlations. Once a choice of composite fields 
has been made, the relative Green's function formalism can be set up, as illustrated 
in Section 2, where a generalized Dyson equation has been derived. On this basis one 
can construct a non-standard perturbation formalism where the basic ingredients are 
the propagators of a subset of the fundamental basis, satisfying linear equations of 
motion [cfr. (3.1)]. By means of the equations of motion and by using the boundary 
conditions related to the definitions of the various Green's functions we have been 
able to derive explicit expressions for these latter [cfr. (3.21), (3.22), (3.45) and 
(3.46)]. However, these expressions can only determine the functional dependence; 
the knowledge of the GF is not fully achieved yet. The reason is that the algebra of 
the field is not canonical. As a consequence, the inhomogeneous terms /^(k) in the 
equations of motion (3.2) and the energy matrix e(k) contain some unknown static 
correlation functions, correlators, that have to be self-consistently calculated. Three 
serious problems arise with the study of the Green's functions: (a) the calculation of 
some parameters expressed as correlation functions of field operators not belonging 
the chosen basis; (b) the appearance of some zero-frequency constants (ZFC) and 
their determination; (c) the problem of fixing the representation where the Green's 
functions are formulated. 

In the Composite Operator Method [26,27] (COM) the three problems (a), (b) 
and (c) are not considered separately but they are all connected in one self-consistent 
scheme. The main idea is that fixing the values of the unknown parameters and of the 
ZFC implies to put some constraints on the representation where the GF are realized. 
As the determination of this representation is not arbitrary, it is clear that there is 
no freedom in fixing these quantities. They must assume values compatible with the 
dynamics and with the right representation. Which is the right representation? 

From the algebra it is possible to derive several relations among the operators. 
We will call algebra constraints (AC) all possible relations among the operators 
dictated by the algebra. This set of relations valid at microscopic level must be 
satisfied also at macroscopic level, when expectations values are considered. Also, we 
note that, in general, the Hamiltonian has some symmetry properties (i.e. rotational 
invariance in coordinate and spin space, phase invariance, gauge invariance,. . . ). 
These symmetries generate a set of relations among the matrix elements: the Ward- 
Takahashi identities [47,48] (WT). 

Now, certainly the right representation must be the one where the relations 
among the operators and the conservation laws are maintained when expectation 
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values are taken; in other words, all the AC and WT are preserved. By imposing 
these conditions we obtain a set of self-consistent equations that will fix the unknown 
correlators, the ZFC and the right representation at the same time. Several equations 
can be written down, according to the different symmetries we want to preserve. A 
large class of self-consistent equations is given by the following equation 

1 1 p + oo 

mW{i)) = ^E^y ^ duC^(Ku) (4.1) 

where the l.h.s. is fixed by the AC, the WT and the boundary conditions compat- 
ible with the phase under investigation and in the r.h.s. the correlation function 
C^t (k, a;) is computed by means of the equation of motion, as illustrated in Sec- 
tion 3. Equations (4.1) generate a set of self-consistent equations which determine 
the unknown parameters (i.e., ZFC and unknown correlators) and, consequently, the 
proper representation [26,27], avoiding the problem of uncontrolled and uncontrol- 
lable decoupling. 



5. Approximation schemes 

The generalized Dyson equation (2.27) is an exact equation and permits, in prin- 
ciple, once the normalization matrix I(i, j) [cfr. Eq. (2.13)], the m-matrix m(i, j) [cfr. 
Eq. (2.14)] and the propagator B(i,j) [cfr. Eq. (2.24)] are known, in the framework 
of the self-consistent scheme outlined in Sections 3 and 4, the calculation of the 
various Green's functions. However, for most of the physical systems of interest the 
calculation of the propagator B(i,j) is a very difficult task and some approximations 
are needed. Various approximate schemes have been proposed. We will summarize 
some of them. 



5.1. The n-pole approximation 

The simplest approximation is based on completely neglecting the dynamical 
part S(k, u) [cfr. (2.28)] of the self-energy. This approximation is largely used in the 
literature [2,6,8-14,19-25,49-62] and is called pole- approximation. In this approx- 
imation we only need the knowledge of the normalization matrix and the m-matrix. 
The constraint (2.9) produces a physics of the solution totally extraneous to the 
complementary physical space, orthogonal to that spanned by the multiplet 
This approximation, or assumption to a larger extent, consists in retaining that one 
can neglect finite life-time effects (i.e., the dynamical part of the self-energy) paying 
attention to the choice of a proper extended operatorial basis, with respect to which 
the self-energy corrections have a small total weight. Indeed, the total weight of the 
corrections is bounded by the thermal average (2.25) involving the residual source 
SJ(i). It is worth noting [63] that the n-pole structure of the various GF corresponds 
to a Dyson-like equation 

<£(k,")= ^ff : (5-1) 
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where the self-energy components S^ 6 (k, uj) have a (n — l)-pole structure. 

5.2. Self-consistent Born approximation 

In order to improve the approximation one needs to take into account self- 
energy corrections by developing some methods to calculate the effects of E(k, uf). 
In the self-consistent Born approximation (SCBA), or non-crossing approximation, 
the many-particle Green's functions, appearing in the expression of S(k, uj) [see 
(2.29)], are calculated by assuming that the fermionic and bosonic modes propagate 
independently. 

By recalling the results given in Section 2, the knowledge of the self-energy 
requires the calculation of the higher-order propagator B Q (i,j) = (Q[5J(i)SJ^(j)}). 
In order to illustrate the approximation, let us consider the case where the basic set 
is of a fermionic type. Then, typically we have to calculate GF of the form 
H R (i,j) = (R[B(i)F(i)F J( (j)B J( (j)]) where F(i) and B(i) are fermionic and bosonic 
field operators, respectively. By means of the spectral representation (3.29) we can 
write 

H R (k,u) = -- / dJ coth^3[# c (k,u/)] (5.2) 

7r uj — uj + ie 2 

where H c (i,j) = (T[B(i)F(i)F^ (j)B^(j)]) is the causal function. In the SCBA we 
approximate H c (i,j) f C (i,j)b c (i,j) where 

f(i,j) = <T[F(i)Ft(j)]> b c ^ 3 ) = (T[B(i)B\j)}) (5.3) 

This approximation has been used in many works (for instance see [64-67]. By 
assuming that the system is ergodic we can use the spectral representations (3.30) 
and (3.54) to obtain 

f(k, W ) = -i <L,'\ l - M ^j + J^-Mf^k, u) ] (5.4) 

7T J_ oc UJ — Uj' + 10 UJ — UJ — 10 



b c (k,uj) = — / duj'[ Jayt ' - , '— \%[b R {k,u)\ (5.5) 

v ' it 1 uj-uj + i5 uj - uj' - i5 i 1 v n v ' 

Use of (5.2)-(5.4) leads to 

H R (k,uj) = - du' — - / d d pd^[f R ( P ,n)} 

kJ-oo uj - uj' + io (27r) d+i J Qb 
%[b R (k - p, J - ft)][tanh ^ + coth (5.6) 

where d is the dimensionality of the system, a is the lattice constant and fie is the 
volume of the Brillouin zone. 

It should be noted that the SCBA can also be applied to the correlation function. 
We start from the expression 

1 f +QO 1 + p~ 0u> 



H R (k,uj) = — du' — -H(k,uj') (5.7) 

2tt / uj-uj' + ie v ' ; K J 
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where H(i,j) = (B(i)F(i)F'(J)B'(j)) is the correlation function. In the SCBA we 
approximate 

H(i,j) = (B{i)F{i)F\j)B\j)) » (F(i)F\j))(B(i)B\j)) (5.8) 
Then, by proceeding in the similar way we arrive to the same expression (5.6). 

5.3. Two-site resolvent approach 

In this subsection we will consider an approximation scheme [68,69], where the 
dynamical part S(k, uj) [cfr. (2.29)] of the self-energy is estimated by a two-site ap- 
proximation in combined use with the resolvent method [70]. In this approximation 
the higher order propagator (2.25) is written as 

B q Ql,w) = F.T.(Q[5J(i)5J\j)\) « B$(u) + aQs)B?(w) (5.9) 

where Bq{oj) is related to level transitions on equal site 

B%(lo) = ±F.T.(R[8J(i,t i )8Ji(i,t j )]) (5.10) 

while Bf(u) is related to transitions across two sites 

B?(u) = ±FT.{R[6J(i,U)5J*(i,t s )]) (5.11) 

The Green's function (2.28) takes the form 

G Q (k,u) = . . - I(k) (5.12) 

v 1 uo - e{k) + t 2 V {uo)a{k) v 1 v 1 

where V(uo) has to be calculated from the definition (2.29) by making use of approx- 
imation (5.9). We give now a brief sketch of the calculation. 

Bq(cu) and B®(u) are computed by expressing 6J(i) in terms of transitions 
among the two-site levels 

6J = J2 a ™®n®rn (5.13) 

nm 

where {$ m } is the complete set of operators for two-site levels. By means of the 
non crossing approximation [70,71], the propagator (Q[$J t (tj)$ m (tj) ( l>J l (t : ,)$ m (t :; )]) 
is expressed in terms of the resolvent R nm (ti — tj) = (Q[$ n (tj)<I > ] ri (t : ,)])i ? where the 
subscript R indicates the reservoir system, i.e., the part Hr of the Hamiltonian 
other than that concerned with two sites. The calculation of the resolvent brings 
to a modification of the original two-site levels by the surroundings. In this scheme 
effects of time delay in the local correlations are treated trough the time-dependent 
modifications of the two-site level transitions, which are included as time-dependent 
local effects in the electron self-energy. 
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6. Conclusions 

In this article we have illustrated an approach to the study of highly correlated 
electronic systems, based on the equation of motion and Green's function method. 
Such an approach is based on two main ideas: (i) propagators of composite opera- 
tors as building blocks at the basis of approximate calculations; (ii) use of algebra 
constrains to fix the representation of the GF in order to maintain the algebraic 
and symmetry properties. This formalism has been applied to the study of several 
models of highly interacting systems, and we refer the interested readers to Ref. [27] 
for an exhaustive bibliography. 

References 

1. P. W. Anderson, Science 235, 1196 (1987). 

2. P. Fulde, Electron Correlations in Molecules and Solids, 3rd ed. ( Springer- Verlag, 
Berlin, 1995). 

3. N. N. Bogoliubov, J. Phys. USSR 11, 23 (1947). 

4. S. M. Dancoff, Phys. Rev. 78, 382 (1950). 

5. R. Zwanzig, in Lectures in Theoretical Physics, edited by W. Britton, B. Downs, and 
J. Downs (Interscience, New York, 1961), Vol. 3, p. 106. 

6. H. Mori, Prog. Theor. Phys. 33, 423 (1965). 

7. H. Umezawa, Advanced Field Theory: Micro, Macro and Thermal Physics (AIP, New 
York, 1993), and references therein. 

8. J. Hubbard, Proc. Roy. Soc. A 276, 238 (1963). 

9. J. Hubbard, Proc. Roy. Soc. A 277, 237 (1964). 

10. J. Hubbard, Proc. Roy. Soc. A 281, 401 (1964). 

11. D. J. Rowe, Rev. Mod. Phys. 40, 153 (1968). 

12. L. M. Roth, Phys. Rev. 184, 451 (1969). 

13. Y. A. Tserkovnikov, Teor. Mat. Fiz. 49, 219 (1981). 

14. Y. A. Tserkovnikov, Teor. Mat. Fiz. 50, 261 (1981). 

15. S. E. Barnes, J. Phys. F 6, 1375 (1976). 

16. P. Coleman, Phys. Rev. B 29, 3035 (1984). 

17. G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57, 1362 (1986). 

18. O. K. Kalashnikov and E. S. Fradkin, Sov. Phys. JETP 28, 317 (1969). 

19. W. Nolting, Z. Phys. 255, 25 (1972). 

20. S. Ishihara, H. Matsumoto, S. Odashima, M. Tachiki, and F. Mancini, Phys. Rev. B 
49, 1350 (1994). 

21. F. Mancini, S. Marra, A. M. Allega, and H. Matsumoto, Physica C 235, 2253 (1994). 

22. F. Mancini, S. Marra, and H. Matsumoto, Physica C 244, 49 (1995). 

23. F. Mancini, S. Marra, and H. Matsumoto, Physica C 250, 184 (1995). 

24. F. Mancini, S. Marra, and H. Matsumoto, Physica C 252, 361 (1995). 

25. A. Avella, F. Mancini, D. Villani, L. Siuraksbina, and V. Y. Yushankhai, 
Int. J. Mod. Phys. B 12, 81 (1998). 

26. F. Mancini and A. Avella, Eur. Phys. J. B 36, 37 (2003). 

27. F. Mancini and A. Avella, Adv. Phys. 53, 537 (2004). 

28. Y. A. Izyumov and B. M. Letfulov, J. Phys.: Condens. Matter 2, 8905 (1990). 



21 



F. Mancini and A. Avella 



29. H. Matsumoto, A. Allega, S. Odashima, and F. Mancini, Physica C 235, 2227 (1994). 

30. F. Mancini, N. B. Perkins, and N. M. Plakida, Phys. Lett. A 284, 286 (2001). 

31. D. Villani, E. Lange, A. Avella, and G. Kotliar, Phys. Rev. Lett. 85, 804 (2000). 

32. A. Avella, R. Hayn, and F. Mancini, Eur. Phys. J. B 37, 465 (2004). 

33. M. Bak and F. Mancini, Physica B 312, 732 (2002). 

34. N. N. Bogoliubov and S. V. Tyablikov, Dokl. Akad. Nauk. USSR 126, 53 (1959). 

35. D. N. Zubarev, Sov. Phys. Uspekhi 3, 320 (1960). 

36. D. N. Zubarev, Non Equilibrium Statistical Thermodynamics (Consultants Bureau, 
New York, 1974). 

37. H. Umezawa, Acta Phys. Hung. 19, 9 (1965). 

38. H. Umezawa, Suppl. of Progr. Theor. Phys. 37-38, 585 (1966). 

39. L. Leplae, F. Mancini, and H. Umezawa, Phys. Rep. 10, 151 (1974). 

40. F. Mancini, Europhys. Lett. 70, 485 (2005). 

41. F. Mancini, Eur. Phys. J. B 45, 497 (2005). 

42. F. Mancini, Eur. Phys. J. B 47, 527 (2005). 

43. A. Avella and F. Mancini, Eur. Phys. J. B 50, 527 (2006). 

44. N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966). 

45. R. Kubo, J. Phys. Soc. Jpn. 12, 570 (1957). 

46. A. Avella, F. Mancini, and E. Plekhanoff, Condens. Matter Phys. 9, ??? (2006). 

47. J. C. Ward, Phys. Rev. 78, 182 (1950). 

48. Y. Takahashi, Nuovo Cimento 6, 370 (1957). 

49. K. W. Becker, W. Brenig, and P. Fulde, Z. Phys. B 81, 165 (1990). 

50. H. Mori, Prog. Theor. Phys. 34, 399 (1965). 

51. A. J. Fedro, Y. Zhou, T. C. Leung, B. N. Harmon, and S. K. Sinha, Phys. Rev. B 46, 
14785 (1992). 

52. N. M. Plakida, V. Y. Yushankhai, and I. V. Stasyuk, Physica C 162-164, 787 (1989). 

53. B. Mehlig, H. Eskes, R. Hayn, and M. B. J. Meinders, Phys. Rev. B 52, 2463 (1995). 

54. J. Beenen and D. M. Edwards, Phys. Rev. B 52, 13636 (1995). 

55. O. K. Kalashnikov and E. S. Fradkin, Phys. Stat. Sol. (b) 59, 9 (1973). 

56. G. Geipel and W. Nolting, Phys. Rev. B 38, 2608 (1988). 

57. W. Nolting and W. Borgiel, Phys. Rev. B 39, 6962 (1989). 

58. A. Lonke, J. Math. Phys. 12, 2422 (1971). 

59. H. Shimahara and S. Takada, J. Phys. Soc. Jpn. 60, 2394 (1991). 

60. P. Kriiger and P. Schuck, Europhys. Lett. 27, 395 (1994). 

61. J. Kondo and K. Yamaji, Prog. Theor. Phys. 47, 807 (1972). 

62. K. Yamaji and J. Kondo, Phys. Lett. A 45, 317 (1973). 

63. F. Mancini, Phys. Lett. A 249, 231 (1998). 

64. N. Plakida and V. Oudovenko, Phys. Rev. B 59, 11949 (1999). 

65. N. M. Plakida, L. Anton, S. Adam, and G. Adam, JETP 97, 331 (2003). 

66. A. Avella, S. Krivenko, F. Mancini, and N. M. Plakida, J. Magn. Magn. Mat. 272, 
456 (2004). 

67. S. Krivenko, A. Avella, F. Mancini, and N. Plakida, Physica B 359-361, 666 (2005). 

68. H. Matsumoto, T. Saikawa, and F. Mancini, Phys. Rev. B 54, 14445 (1996). 

69. H. Matsumoto and F. Mancini, Phys. Rev. B 55, 2095 (1997). 

70. Y. Kuramoto, Z. Phys. B 53, 37 (1983). 

71. N. Grewe, Z. Phys. B 53, 271 (1983). 



22 



